function S = gmm_gibbs_demo(K, n)
% A simple demo on fmm_gibbs
%
%   S = gmm_gibbs_demo(K, n);
%
%       This function demonstrates the Gibbs sampling over finite mixture
%       model with 2D Gaussian component models.
%
%       Inputs:
%       - K:    the number of components
%       - n:    the number of samples from each component.
%
%       Outputs:
%       - S:    the consoliated samples
%

% Created by Dahua Lin, on Nov 16, 2010
%

%% prepare model

disp('preparing model ...');

d = 2;

prior = gaussd.from_mp(0, udmat(d, 6^2), 'ip');
mu = prior.sample(K);
sigma = udmat(d, 1);

gm = gaussgm_gp(prior, [], sigma);



%% synthesize data

disp('synthesizing data ...');

X = zeros(d, n*K);
for k = 1 : K
    X(:, (k-1)*n + (1:n)) = gm.sample(mu(:,k), n);
end

%% do sampling

disp('Doing sampling ...');

R = fmm_gibbs(gm, X, K, 500, 'BurnIn', 1, 'Interval', 2);


%% consoliate samples

disp('Consolidating samples ...');
S = consolidate_mmsample(R, 1e-2 * (1/K));

for k = 1 : numel(S.Id)    
    theta = S.Theta(:, k);
    fprintf('[%d]: (w = %.4g): [%s]\n', ...
        S.Id(k), S.W(k), num2str(theta', ' %.4f'));    
end

disp(' ');
G = gaussd.from_mp(S.Theta, sigma);


%% visualize

% plot data

figure;
plot(X(1,:), X(2,:), 'b.', 'MarkerSize', 3);

hold on;
plot_ellipse(G, 1, 'g-', 'LineWidth', 2);
plot_ellipse(G, 3, 'g-');

axis equal;

% plot sample traces

Thetas = cat(3, R.samples.Theta);

for k = 1 : K    
    T = squeeze(Thetas(:, k, :));
    hold on;
    plot(T(1,:), T(2,:), 'r-');    
    hold on;
    plot(S.Theta(1,k), S.Theta(2,k), 'm+', 'MarkerSize', 15);
end


